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ABSTRACT 


A  finite  element  approach  is  used  to  obtain  equations  for  large 
displacement  and  stability  analysis  via  the  principle  of  virtual  work. 

The  equations  are  then  written  in  incremental  form.  It  is  found  that  the 
inclusion  of  the  nonlinear  strain  displacement  terms  results  in  two  addi¬ 
tional  ’initial  stiffness’  matrices.  The  first  is  the  well  known  'initial 
stress'  matrix.  An  alternate  general  method  is  proposed  for  the  derivation 
of  the  'initial  stress'  matrix.  The  second  matrix  is  a  function  of  previous 
displacements  and  is  referred  to  here  as  the  'initial  displacement'  matrix. 
It  is  obtained  as  a  result  of  writing  the  nonlinear  strain  displacement 
relations  in  incremental  form.  The  'initial  displacement'  matrix  was 
found  to  be  of  the  same  order  as  the  'initial  stress'  matrix  but  appears 
not  to  have  been  previously  recognized  in  finite  element  analysis. 

Results  obtained  for  a  simple  truss  and  an  arch  problem  agreed  with 
results  in  the  literature.  Analysis  without  the  'initial  displacement' 
matrix  resulted  in  an  overestimate  of  the  buckling  loads  by  factors  of 
1.8  and  2  for  the  truss  and  arch  respectively.  A  linear  eigenvalue  analysis 
resulted  in  an  overestimate  of  the  buckling  load  by  factors  of  5  and  2 
respectively. 

It  was  concluded  that  'initial  displacement'  stiffness  matrices 
should  be  included  in  the  finite  element  analysis  of  large  deformation 
and  instability  problems. 


iii 


NOTATION  AND  DEFINITIONS 


{a} 

[B] 

CD] 

Ce> 

[f(u,x>] 

[k(o)] 

[k(1)] 

{P> 

[T(u,x)] 

u,v 

{u} 


V 

{x} 

9 

A 

{a} 

(u,x) 


generalized  coordinates,  coefficients  of  displacement 
polynomials  =  [a]{u>. 

differential  operator  that  transforms  generalized 
coordinates  to  increments  of  strains.  Ae^  =  [B]A{a>. 
constitutive  relation,  strain  to  stress  transformation  matrix, 
strain  vector. 

displacement  strain  relation. 

stiffness  matrix  due  to  nonlinear  strain-displacement  relation, 
stiffness  matrix  due  to  initial  stress, 
load  vector. 

local  to  global  coordinate  transformation  matrix, 
axial  and  normal  displacements,  respectively, 
displacement  vector,  referred  to  local  coordinate  system 
when  not  preceded  by  [T(u,x)]. 
volume . 

position  vector, 
a  virtual  change, 
incremental  change, 
stress  vector. 

function  of  displacement  and  geometry. 


Subscripts 

i 

o,l 


subscript  denoting  matrix  row  or  column, 
subscript  denoting  initial  and  final  geometry. 
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INTRODUCTION 

Large  displacement  analysis  by  the  finite  element  method  was  first 

proposed  by  Turner  et  al.^~  'Initial  stress'  matrices  were  suggested  to 

account  for  the  effect  of  initial  stress  in  truss  and  plane  stress  assemblies. 

Subsequent  work  on  the  derivation  of  the  'initial  stress'  matrices  for  other 
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elements  were  reported  by  Argyris  et  al,  Gallagher  et  al  and  Kapur  and 
4  5 

Hartz.  In  a  recent  work,  Martin  obtained  the  'initial  stress'  matrix  by 
considering  the  potential  energy  of  the  element.  Nonlinear  displacement-strain 
relations  from  the  classical  theory  of  elasticity  were  used. 

This  paper  examines  large  displacement  element  behavior  by  means  of  the 
principle  of  virtual  work.  An  alternate  general  method  is  obtained  for  finding 
the  'initial  stress'  matrix.  All  previous  authors  have  neglected  the  contribu¬ 
tion  of  initial  displacements  to  the  element  stiffness  matrices.  This  paper 
considers  the  effect  of  including  the  initial  displacement  terms  and  shows 
that  it  results  in  'initial  displacement'  matrices  which  are  of  the  same  order 
as  the  'initial  stress'  matrices. 

Numerical  results  are  then  presented  to  justify  the  introduction  of  the 
new  'initial  displacement'  stiffness  matrix. 


PRINCIPLE  OF  VIRTUAL  WORK 

The  equations  for  large  deformation  analysis  are  derived  by  the 
principle  of  virtual  work.  The  general  equations  are  considered  without 
regard  to  subscription  for  elements  and  nodal  points.  The  subdivision  and 
expression  of  the  formulation  into  finite  elements  and  assembly  into  large 
master  stiffness  and  other  matrices  follow  the  usual  procedure. 

In  general  the  strain  in  the  element  may  be  described  as  a  nonlinear 
function  of  the  current  geometry,  previous  displacements  and  the  displacement 
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at  the  nodes  {u>. 


{e}  =  f(u,x)  (1) 

The  strain  due  to  a  virtual  displacement  3{u}  is 

3{e}  =  [  |f  (u,x)][T(u,x)]3{u)  (2) 

where  the  transformation  matrix  [T(u,x)]  has  been  included  to  indicate 
transformation  from  local  to  global  coordinate  systems.  The  brackets  (u,x) 
are  used  to  indicate  the  dependence  of  the  matrix  on  current  nodal  displace¬ 
ments  and  geometry  respectively. 

By  the  principle  of  virtual  work 

|  3LuJ[T(u,x)]T  [  |f  (u,x)]T{o>dV  =  3LuJ{P}  (3) 

V 


so  that 

|  [T(u,x)]T  [  |f  (u,x)]T{o}dV  =  {P>  (4) 

V 

where  {P}  is  the  load  vector  at  the  nodes  and  should  be  understood  to  include 
the  effect  of  uniform  loading. 

This  is  the  equation  describing  the  finite  deformation  behavior  of  a 
solid.  The  dependence  of  the  various  matrices  on  geometry  and  previous  dis¬ 
placements  are  to  be  particularly  noted. 


ITERATIVE  SOLUTIONS 

6 

Equation  (4)  has  been  solved  by  Bogner  e't  al  by  the  use  of  energy  search 
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procedures.  Turner  et  al  and  subsequent  work  summarized  by  Zienkiewicz 

however  have  sought  to  linearize  Eq.  (4)  and  account  for  the  large  displace¬ 
ment  by  iterative  solutions  for  the  final  geometry.  A  symmetric  stiffness 
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element  based  on  final  geometry  was  used.  It  can  be  seen  that,  for  an 
elastic  structure,  this  is  equivalent  to  replacing  the  stress  vector  by 

{0}  =  [E]  [  ||  (u,x)]1CT(u,x)]1{u}  (5) 

where  subscript  1  denotes  values  at  the  end  of  the  loading.  It  is  the 
author's  opinion  that  a  better  approximation  would  result  from 

{0}  =  |  [E]  ([  ||  (o,x)]o[T(o,x)]o  +  [  ||  (u,x)]1[T(u,x)]1){u}  (6) 

where  subscript  o  denotes  initial  values. 

Thus  the  linearized  iterative  solutions  may  be  interpreted  as  application 
of  integration  methods  where  strain  rates  at  the  beginning  and  end  of  the 
integration  are  used.  Equation  (6)  is  a  higher  order  formula  and  still  more 
refined  methods  may  be  introduced  which  in  the  limit  will  tend  to  the 
incremental  solution. 

INCREMENTAL  SOLUTIONS 

Equation  (4)  is  a  nonlinear  equation  and  must  be  solved  by  iteration. 

It  is  often  more  advantageous  to  solve  a  series  of  linear  equations.  This 
has  given  rise  to  the  incremental  piece-wise  linear  solutions. 

To  obtain  the  appropriate  equations,  Eq.  (4)  is  subjected  to*  a  small 
perturbation  and  to  first  order 

|  A[T(u,x)]T  [  ||  (u,x)]T{o}dV  +  J  [T(u,x)]TA[  ||  (u,x)]T{0}dV 
V  V 

+  J  [T(u,x)]T  [  ||  (u,x)]T  [D]  [  ||  (u,x)]  [T(u,x)]A{u}dV  =  A{P}  (7) 
V 

where  [D]  is  the  linearized  stress-strain  relation  for  the  increment. 
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The  first  term  in  Eq.  (7)  gives  the  so-called  'geometric  stiffness' 

2  5 

matrix,  the  second  term  will  be  shown  to  be  the  'initial  stress'  matrix 
and  the  third  term  is  the  element  stiffness  matrix  which  also  contains  the 
'initial  displacement'  matrix.  Following  previous  works,  *  the  geometric 
stiffness  is  neglected  in  the  present  study. 

INITIAL  STRESS  MATRIX 

We  now  show  a  method  for  obtaining  the  'initial  stress'  matrix  from  the 
second  term  in  Eq.  (7).  In  this  approach,  separate  matrices  are  first  formed 
for  each  row  of  the  initial  stress  vector  {a}. 

3f 

Hence  the  matrix  term  due  to  the  ith  column  of  [^j]  and  ith  row  of 
{o}  is  given  by 


(9) 


(10) 


(ID 
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3{B}T 


|  a.[T]T  A{|^}  dV  =  |  a.[T]T  [a]T  C^p]  [a][T]A{u}dV 

ir  ^  ii 


=  [k^1)]A{u) 


(12) 


The  full  'initial  stress'  matrix  may  then  be  obtained  by  summing  all  the 
terms  for  all  the  stresses.  It  is  noted  that  the  'initial  stress'  matrix  is 
dependent  on  the  current  stress  and  the  current  geometry.  This  procedure 
yields  the  same  ' initial  stress '  matrix  as  the  procedures  suggested  in 
References  3-5.  It  is  described  here  for  its  compactness  and  generality. 

For  example,  it  may  be  used  for  any  element  with  any  degree  of  nonlinearity 
included  in  the  strain  displacement  equations. 


INITIAL  DISPLACEMENT  MATRIX 


We  consider  the  effect  of  initial  displacements  on  the  third  term  of 
Eq.  (7).  A  second  order  strain  displacement  relation  (usually  of  rotation) 
results  in  a  first  order  current  displacement  term  in  the  matrix  [jjj  (u,x)]. 
If  [D]  is  essentially  constant  throughout  the  analysis,  pre-multiplication 
by  [D]  results  in  terms  of  the  order  of  the  total  stress.  Hence  the  third 
term  in  Eq.  (7)  yields  additional  terms  which  are  of  the  same  order  as  the 
'initial  stress'  matrix.  These  terms  in  the  element  stiffness  are  referred 
to  here  as  the  'initial  displacement'  matrix  in  analogy  with  the  'initial 
stress'  matrix.  As  far  as  the  author  is  aware,  the  'initial  displacement' 
matrix  has  been  ignored  in  all  previous  matrix  finite  element  analysis. 

Its  effect  is,  of  course,  implicitly  included  in  the  analysis  of 
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The  initial  stress  matrix  and  the  element  stiffness  matrix  is  derived  for 
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the  beam-column.  The  notation  of  Reference  5  will  be  adopted  here  with  the 
same  assumption  of  linear  and  cubic  displacements  in  the  axial  and  normal 
displacements  u  and  v. 


e 


du  1^  fdv-) 2 

dx  2  W 


u  =  a  +  a.x 

o  1 

v  =  b  +  b,x  +  b0x2  +  b0x^ 

o  1  2  3 


(13) 


(14) 


separating  into  membrane  and  bending  components. 
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A{e}  = 


0  1  o  p-  2x  p-  3x2  p- 

dx  dx  dx 


-2y 


-6xy 


{a} 


(17) 


=  [B] A  {a}  =  [B]  [a]  A  {u} 

Now  only  considering  the  membrane  strain  effect,  i.e.,  only  taking  the  first 
column  of  matrix  [B] , 


3{B), 


m 


9{a) 


■  < 


0 

N 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

dv 

dx 

>  - 

0 

0 

0 

1 

2x 

3x2 

2x  — 
dx 

0 

0 

0 

2x 

4x2 

6x3 

3x2  p- 
dx 

\ 

J 

0 

0 

0 

3x2 

6x3 

9x4 

3{a) 


By  Eq.  (12)  the  initial  stress  matrix  becomes 


(18) 


Cki1>] 


=  I 


o  [T]TCa]T 
m 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 


0 

0 

0 

1 


0  2x 


0 

0 

0 

2x 

4x2 


0 

0 

0 

3x2 
6x  3 


0  3x2  6x3  9x4 


[a][T]dV 


(19) 


This  is  identical  with  the  expression  derived  by  Martin.' 


(o)- 


The  third  term  gives  rise  to  the  element  stiffness  matrix  [kv  ']  for 


elastic  behavior  and  neglecting  second  order  terms,  we  have  the  element 


8 


stiffness  matrix 


0 

0  1 


DO 


’I 


[T]  [a]1  E 


0 

0 

0 

0 


dv 

dx 


0 

0  0 


2x 


3x 


dv 

dx 

gdv 

dx 


0  0  4y2 

0  0  12xy2  36xy2 


Ca][T]dV 


(20) 


There  are  linear  coefficients  in  displacements  in  the  stiffness  matrix  [k^°*]. 
These  linear  coefficients  form  the  'initial  displacement'  matrix. 


INSTABILITY 

The  instability  that  is  sought  exists  as  stationary  points  in  the  solution 
of  Eq.  (4).  It  is  only  when  the  incremental  equation  (7)  is  used  that  we 
obtain  the  eigenvalue  problem. 

From  this  standpoint,  we  discern  several  levels  of  approximation  in  the 
5  2 

literature.  Martin  and  Argyris  suggested  a  straightforward  use  of  the 
'initial  stress'  matrix  by  scaling  as  in  the  classical  eigenvalue  approach. 

3 

Gallagher  et  al  recognized  the  need  for  the  correct  evaluation  of  the  stress 
terms  in  the  'initial  stress'  matrix.  This  was  achieved  by  iteration.  An 
improved  level  of  approximation  may  be  achieved  by  inclusion  of  the  'initial 
displacement'  matrix  in  an  incremental  large  deformation  analysis.  Unless 
the  mode  of  instability  is  known,  the  master  stiffness  equations  must  be 
examined  for  instability  after  each  increment  of  load.  That  is  to  say 
unsymmetric  modes  may  not  be  obtained  from  a  symmetric  pre-buckling 
displacement  pattern  in  an  incremental  solution. 
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COMPUTER  PROGRAM 

A  general  purpose  finite  element  computer  program  was  modified  to  allow 
a  power  sweep  for  the  lowest  eigenvalue  and  eigenvector  after  each  increment 
of  load.  The  geometry  was  updated  after  each  increment  of  load.  The  analysis 
was  load  controlled  as  opposed  to  deflection  controlled.  It  was  thought  that 
such  an  approach  would  allow  a  search  for  asymmetric  modes  in  cases  with  non¬ 
proportional  loading.  For  this  reason,  no  results  were  obtained  for  the 
post-buckled  region.  Stiffness  matrices  were  developed  for  the  truss  and 
beam-column  element.  The  arch  was  represented  by  a  series  of  straight 
beam-column  elements. 

NUMERICAL  EXAMPLES 

Results  were  obtained  for  a  simple  truss  and  an  arch  example.  Figure  1 

0 

shows  the  single-bar  structure  analyzed  by  Mallett  and  Berke.  The  results 
from  the  present  theory  agree  with  that  of  Reference  8.  The  effect  of  neglect¬ 
ing  the  'initial  displacement'  matrix  is  also  shown.  The  linear  instability 
load  is  included  to  show  the  effect  of  neglecting  changes  in  geometry  and 
stress. 

9 

Figure  2  shows  an  arch  tested  by  Gjelsvik  and  Bodner.  The  experimental 
and  predicted  buckling  loads  agree  when  the  analysis  is  performed  with  16 
elements.  Again,  results  obtained  without  the  'initial  displacement'  matrix 
demonstrate  the  need  for  this  matrix. 

DISCUSSION 

The  results  for  the  arch  show  an  unexpected  variation  with  the  number  of 
elements.  Two  factors  contribute  to  this.  First  of  all  the  nature  of  the 
point  load  can  be  expected  to  cause  rapid  changes  of  displacements  near  its 
point  of  application.  The  second  factor  is  that  the  curved  arch  is  being 
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approximated  by  a  straight  beam  column. 

Calculations  with  32  elements  showed  yet  another  effect  of  using  the 
beam  column  approximation.  This  was  that  the  equations  of  equilibrium  showed 
poor  convergence.  The  cause  of  this  may  best  be  explained  by  noting  that  the 
representation  of  the  correct  arch  effect  by  a  beam  element  depends  on  the 
axial  component  of  the  relative  vertical  displacements  between  the  nodes. 

For  the  element  adjacent  to  the  mid-point  of  the  arch,  both  the  slope  of  the 
element  and  the  relative  vertical  displacements  are  obtained  as  a  result  of 
small  differences  of  large  quantities.  Hence,  appreciable  round-off  errors 
can  be  expected  in  solutions  with  small  elements.  This  was  the  reason  for 
abandoning  the  analysis  with  32  elements. 

The  choice  of  a  linear  axial  displacement  which  resulted  in  a  constant 
axial  stress  was  not  a  good  one.  The  second  order  strain  terms  resulted  in 
a  saw-toothed  axial  stress  distribution  along  the  arch.  For  the  present,  the 
saw-toothed  stress  distribution  was  removed  by  assuming  that  the  stress  state 
throughout  each  arch  element  was  given  by  its  centroidal  value.  Though  the 
results  are  no  longer  'consistent*  within  the  principle  of  virtual  work,  this 
stop-gap  element  allows  an  investigation  of  the  'initial  displacement'  matrix. 

Reference  8  only  required  4  elements  to  obtain  good  results.  It  is 
thought  that  this  is  due  to  the  element  used  there  which  always  satisfied 
axial  equilibrium.  It  may  be  that,  with  the  use  of  an  improved  axial  displace¬ 
ment  function,  the  same  results  could  be  achieved  here  with  fewer  elements. 

The  results  in  the  above  examples  demonstrate  the  need  for  including  the 
'initial  displacement'  matrix  in  a  large  deformation  analysis.  The  structures 
tested,  however,  are  noted  for  their  sensitivity  to  changes  in  geometry.  It 
may  well  be  that,  in  general,  initial  displacements  do  not  play  such  a  large 


role. 
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NOTE  ON  LARGE  DISPLACEMENT  ANALYSIS  BY  ASSUMED  MODES 

The  'initial  displacement'  matrix  has  been  obtained  here  and  shown  to  be 
as  important  as  the  'initial  stress'  matrix.  This  appears  to  be  a  new  effect 
in  the  context  of  finite  element  analysis.  However,  it  is  important  to  note 
that,  within  the  larger  framework  of  the  literature  on  large  displacement  anal¬ 
ysis  in  continuum  mechanics,  the  role  of  initial  displacements  is  not  new  and 
has  always  been  recognized.  It  is  particularly  interesting  to  regard  the 
finite  element  method  as  a  special  case  of  the  Ritz  assumed  modes  technique 
in  which  the  modes  are  written  in  terms  of  nodal  point  displacements.  A 
large  body  of  literature  already  exists  on  the  use  of  the  Ritz  assumed  modes 
technique  for  the  solution  of  large  displacement  and  instability  problems. 

For  example,  in  a  summary  of  Koiter's  work^  dealing  with  post-buckling  be¬ 
havior,  Budiansky"^  has  presented  equations  which  can  be  seen  to  include 
Eqs.  (4)  and  (7)  as  special  cases  (from  Eqs.  (3)  and  (4)  in  Reference  11). 

In  general,  the  equations'^  do  not  lead  to  such  a  simple  and  clear-cut  matrix 
as  the  'initial  displacement*  matrix  found  here.  The  difference  between  the 
present  work  and  References  10  and  11  is  that  as  a  first  approximation,  the 
latter  works  neglect  the  effect  of  displacements  prior  to  instability.  This 
displacement  has  been  shown  to  be  important  in  the  examples  considered  here. 

CONCLUSIONS  AND  RECOMMENDATIONS 

1.  A  compact  and  general  method  has  been  suggested  for  finding  the  'initial 
stress'  matrix  for  finite  element  analysis. 

2.  Initial  displacements  were  shown  to  contribute  an  'initial  displacement* 
matrix  to  the  stiffness  matrix  of  an  element  in  a  large  deformation 
analysis.  This  matrix  was  seen  to  be  of  the  same  order  as  that  of  the 
'initial  stress'  matrix. 
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3.  Numerical  examples  were  presented  which  justify  the  introduction  of 
the  initial  displacement  matrix. 

4.  The  effect  of  including  the  'initial  displacement'  matrix  should  be 
investigated  over  a  wider  variety  of  elements  and  material  behavior. 
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